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ABSTRACT 

The dynamical features of the irregular satellites of the giant planets argue against an 
in-situ formation and are strongly suggestive of a capture origin. Since the last detailed 
investigations of their dynamics, the total number of satellites have doubled, increasing 
from 50 to 109, and almost tripled in the case of Saturn system. We have performed 
a new dynamical exploration of Saturn system to test whether the larger sample of 
bodies could improve our understanding of which dynamical features are primordial 
and which are the outcome of the secular evolution of the system. We have performed 
detailed N-Body simulations using the best orbital data available and analysed the 
frequencies of motion to search for resonances and other possible perturbing effects. 
We took advantage of the Hierarchical Jacobian Symplectic algorithm to include in 
the dynamical model of the system also the gravitational effects of the two outermost 
massive satellites, Titan and Iapetus. Our results suggest that Saturn's irregular satel- 
lites have been significantly altered and shaped by the gravitational perturbations of 
Jupiter, Titan, Iapetus and the Sun and by the collisional sweeping effect of Phoebe. 
In particular, the effects on the dynamical evolution of the system of the two massive 
satellites appear to be non-negligible. Jupiter perturbs the satellites through its direct 
gravitational pull and, indirectly, via the effects of the Great Inequality, i.e. its almost 
resonance with Saturn. Finally, by using the Hierarchical Clustering Method we found 
hints to the existence of collisional families and compared them with the available 
observational data. 

Key words: planets and satellites: formation, Saturn, irregular satellites - methods: 
numerical, N-Body simulations - celestial mechanics 



1 INTRODUCTION 

The outer Solar System is inhabited by different minor bod- 
ies populations: comets, Centaurs, trans-neptunian objects 
(TNO), Trojans and irregular satellites of the giant planets. 
Apart from comets, whose existence had long been known, 
the other populations have been discovered in the last cen- 
tury and extensively studied since then. Our comprehension 
of their dynamical and physical histories has greatly im- 



proved (see Jewitt (20081; Morbidelli ( 2008 I for a general 



review on the subject) and recent models of the formation 
and evolution of the Solar System seem to succeed in ex- 



plaining their evolution and overall structure (see Gomes et 
ah] p005l); [Morbidelli etld~1([2005|);[Tsiganis et aT|||2005} for 



details). There is still no clear consensus, however, on the ori- 



gin of the irregular satellites. It's widely accepted that their 
orbital features are not compatible with an in situ forma- 
tion, leading to the conclusion they must be captured bodies. 
At present, however, no single capture model is universally 



accepted (see Sheppard ( 2006 1 and Jewitt & Haghighipour 
[] pOOTt for a general review). It is historically believed that 
the satellite capture occurred prior to the dissipation of the 
Solar Nebula ( |Pollack et al.|1979[ ), since the gaseous drag is 
essential for the energy loss process that leads to the cap- 
ture of a body in a satellite orbit while crossing the sphere 
of influence of a giant planet. 

A direct consequence and implicit assumption of this model 
is that no major removal of irregular satellites took place 
after the dissipation of the nebular gas for the actual satel- 
lite systems to be representative of the gas-captured ones. If 
instead this was the case, we are left with two possibilities: 
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the capture efficiency of the gas-drag mechanism was by 



2 Turrini et al. 



Table 1. Mean orbital elements for the 35 irregular satellites of Saturn computed with Model 1. For each orbital element, we report the 
minimum and maximum values attained during the simulations. 



Mean Max Min Mean Max Min Mean Max Min 

Satellite a a a e e e i i i 
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far higher than previously assumed in order to have enough 
satellites surviving to the present time in spite of the re- 
moval, or 

• a second phase of capture events, based on different 
physical mechanisms, took place after the gas dissipation. 

The same issues apply also to the original formulation of the 
so called Pull-Down scenario (Hcppcnh eimer fc Porco|1977[ ), 
also based on the presence of the nebular gas but locating 
the capture events during the phase of rapid gas accretion 
and mass growth of the giant planets. 

Recently, the plausibility of gas-based scenarios has been put 
on jeopardy. The comparative study performed by |Jewitt fc| 
|Sheppard~| ( |2005[ > pointed out that the giant planets possess 
similar abundances of irregular satellites, once their appar- 
ent magnitudes are scaled and corrected to match the same 
geocentric distance. Due to the different formation histories 
of gas and ice giants, the gas-based scenarios cannot supply 



a convincing explanation to this fact. 



The Nice model (Gomes et al. 2005 Morbidelli et al. 2005 



Tsiganis et al. 2005) instead undermines the physical rele- 



vance of the gas-based scenarios. Formulated to explain the 
present orbital structure of the outer Solar System, it postu- 
lates that the giant planets formed (or migrated due to the 
interaction with the nebular gas) in a more compact con- 
figuration than the current one. Successively, due to their 
mutual gravitational interactions, the giant planets evolved 
through a phase of dynamical rearrangement in which the 
ice giants and Saturn migrated outward while Jupiter moved 
slightly inward. The dynamical evolution of the giant planets 
was stabilised, during this phase, by the interaction with a 
disc of residual planetesimals. During the migration process, 
a fraction of the planetesimals can be captured as Trojans 
( Bottke et aL|20 08) while the surviving outer planetesimal 
disk would slowly settle down as the present Kuiper Belt. 
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Table 2. Mean orbital elements for the 35 irregular satellites of Saturn computed with Model 2. For each orbital element, we report the 
minimum and maximum values attained during the simulations. 



Mean Max Min Mean Max Min Mean Max Min 
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The perturbations of the migrating planets on the planetes- 
imals would have also caused the onset of the Late Heavy 
Bombardment on the inner planets. The detailed descrip- 
tion of this model and its implications are given in Gomes| 



et al. (12005 ); Morbidelli et al. (20051; Tsiganis et al 



(2005) 



and in |Morbidelli et ah (20071. A major issue with the orig- 
inal formulation of the Nice model was the ad hoc choice 
of the orbits of the giant planets after the dispersal of the 



nebular gas. Work is ongoing to solve the issue (Morbidelli 



et al. 20071 by assuming that migration of planets by in- 



teraction with the Solar Nebula could lead to a planetary 
resonant configuration subsequently destroyed by the inter- 
action with the planetesimal disk. 

The Nice model has two major consequences concerning the 
origin of the irregular satellites. Due to the large distance 
from their parent planets and their consequently looser grav- 
itational bounds, no irregular satellites previously captured 



could have survived the phase of violent rearrangement tak- 
ing place in the Solar System ( ] Tsiganis et al.|2"505| ). In ad- 
dition, during the rearrangement phase planetesimals where 
crossing the region of the giant planets with an intensity 
orders of magnitude superior to the one we observe in the 
present Solar System (Tsiganis et al. 2005), favouring cap- 
ture by different dynamical processes. 

In other words, the Nice model doesn't rule out gas-based 
mechanisms for irregular satellites capture but it implies 
that they could not have survived till now due to the vio- 
lent dynamical evolution of the planets. The only exception 
to the former statement could be represented by the Jo- 
vian system of irregular satellites, since Jupiter had a some- 
how quieter dynamical evolution and its satellites were more 
strongly tied to the planet due to its intense gravity. For the 
other giant planets, however, the chaotic orbital evolution 
created a favourable environment for the capture of new ir- 
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regular satellites by three body processes (i.e. gravitational 
interactions during close encounters and collisions). These 
trapping processes don't significantly depend on the choice 
of the initial orbital configuration of the planets in the Nice 
model: they rely mainly on the presence of a residual disk of 
planetesimals which both stabilises the dynamical evolution 
of the giant planets and supplies the bodies that can become 
irregular satellites. 

In this paper we concentrate on the study of the dynam- 
ical features of Saturn's irregular satellites, stimulated by 
the unprecedented data gathered by the Cassini mission on 
Phoebe. The overall structure of the satellite system can in 
fact provide significant clues on its origin and capture mech- 
anism. This dynamical exploration will be the first step to 
test the feasibility and implications of the collisional cap- 
ture scenario which was originally proposed by |Colombo fc| 
Franklin ( 1971 1, which will be the subject of a forthcoming 



paper. The work we will present is organised as follows: 

• in section [2] we will describe the setup we used to simu- 
late the dynamical history of the irregular satellites of Sat- 
urn through N-Body integrations and present the proper 
orbital elements we computed; 

• in section [3] we will analyse the dynamical evolution of 
the satellites and search for resonant or chaotic behaviours; 

• in section[4]we will present the updated collisional prob- 
abilities between pairs of irregular satellites. We will also 
estimate the impact probability of bodies orbiting close to 
Phoebe to explore the relevance of post-capture impacts in 
sculpting the system; 

• in section [5] we will apply the Hierarchical Clustering 
Method described in Zappala et al. ( 1990 1994 I to search for 



dynamical families and compare the results with the avail- 
able observational data; 

• finally, in section [6] we will put together all our results 
to draw a global picture of Saturn's system of irregular satel- 
lites. 



2 COMPUTATION OF THE MEAN ORBITAL 
ELEMENTS 

Our study of the dynamics of the irregular satellites start 
from the computations of approximate proper orbital ele- 
ments. Proper elements may be derived analytically from 



the nonlinear theory developed by Milani & Knezevic ( 1994 1 



or they can be approximated by the mean orbital elements, 
which can be numerically computed by averaging the oscu- 
lating orbital elements over a reasonably long time interval 
which depends on the evolution of the system. 
To compute the mean orbital elements of Saturn's irregular 
satellites we adopted the numerical approach and integrated 
over 10 8 years the evolution of an N-body dynamical system 
composed of the Sun, the four giant planets and the satel- 
lites of Saturn in orbit around the planet. For the satellites, 
we adopted two different dynamical setups: 

(i) Model 1: 35 irregular satellites treated as test particles 

(ii) Model 2: 35 irregular satellites treated as test parti- 
cles, in addition Titan and Iapetus treated as massive bodies 

Model 1 basically reproduces the orbital scheme used by 
other authors in pr evious works (|Carruba et ah] |2002| 
Nesvorny et al.|2003| |Cuk fc Burns|2004| >. Model 2 includes 
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Figure 1. Comparison in the a— i plane between the mean orbital 
elements of Saturn's irregular satellites computed with Model 2 
using HJS and RADAU algorithms (upper plot), Model 1 and 
Model 2 using HJS algorithm (middle plot) and Model 1 with 
standard and strict double precision (lower plot). The vertical and 
horizontal bars show the variation ranges of the elements in the 
simulations. Distances are expressed in 10 6 km while angles are 
expressed in degrees. In the simulation with RADAU algorithm 
we considered only the 26 irregular satellites known till 2005. 



the perturbing effects of the two outermost regular satellites 
of the Saturn system. We considered this alternative model 
to evaluate their contribution to the global stability and to 
the secular evolution of the irregular satellites. 
In order to reproduce accurately the orbits of the satellites, 
the integration timestep has to be less than the shortest 
period associated to the frequencies of motion of the sys- 
tem: the common prescription for computational celestial 
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Figure 2. Comparison in the a— e plane between the mean orbital 
elements of Saturn's irregular satellites computed with Model 2 
using HJS and RADAU algorithms (upper plot), Model 1 and 
Model 2 using HJS algorithm (middle plot) and Model 1 with 
standard and strict double precision (lower plot). The vertical and 
horizontal bars show the variation ranges of the elements in the 
simulations. Distances are expressed in 10 6 km while angles are 
expressed in degrees. In the simulation with RADAU algorithm 
we considered only the 26 irregular satellites known till 2005. 



mechanics is to use integration timesteps about 1/20 of the 
fastest orbital period, usually that of the innermost body. 
The timesteps we used in our simulations with Models 1 and 
2 were respectively: 

• 20 days (1/22 of the orbital period of Kiviuq) 

• 0.75 day (1/20 of the orbital period of Titan) 



Figure 3. Comparison between the mean orbital elements of Sat- 
urn's irregular satellites computed with Model 2 using the HJS 
algorithm and those from the JPL Solar System Dynamics web- 
site. The comparison is shown in the a — i (upper graph) and a — e 
(lower graph) planes. The vertical and horizontal bars show the 
variation ranges of the elements in our simulation. Distances are 
expressed in 10 6 km while angles are expressed in degrees. 



The timestep used for Model 2 greatly limited the length of 
the simulations by imposing a heavier computational load. A 
time interval of 10 8 years was a reasonable compromise be- 
tween the need of accuracy in the computation of proper ele- 
ments and CPU (i.e. computational time) requirements: this 

((20031. 



time interval has been used also by Nesvorny et al 



We used it for both Models in order to be able to compare 
their results. 

The initial osculating orbital elements of the irregular satel- 
lites and of the major bodies have been derived respectively 
from two different ephemeris services: 

• Natural Satellites of the I A U Minor Planet Cente^\ 

• Horizons of the NASA Jet Propulsion Laboratorjj^] 

The reference plane in all the simulations was the J2000 
ecliptic and the initial elements of all bodies referred to the 
epoch 30 January 2005. 

The numerical simulations have been performed with the 
public implementation of the HJS (Hierarchical Jacobi Sym- 
plectic) algorithm described in |Beust| ( |2003| ) and based on 



http : //cf a- www .harvard. edu/ iau/NatSats 



http : //ssd. jpl .nasa. gov/?horizons 
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the SWIFT code (jLevison fc Duncanp.994). The symplec 



tic mapping scheme of HJS is applicable to any hierarchi- 
cal system without any a priori restriction on the orbital 
structure. The public implementation does not include the 
effects of planetary oblateness or tidal forces which are not 
considered in our model^] The HJS algorithm proved itself 
a valuable tool to study the dynamical evolution of satel- 
lite systems since its symplectic scheme supports multiple 
orbital centres. It easily allows to include in the simulations 
both Titan and Iapetus which were not taken into account 
in the previous works on irregular satellites. 
To test the reliability of HJS, we performed an additional 



run based on Model 2 using the RADAU algorithm (Ever- 
hart||1985[ > incor porated in the M ERCURY 6.2 package by 



J.E. Chambers (Chambers 19991. Due to the high compu- 



tational load of RADAU algorithm compared to symplectic 
mapping, we limited this simulation to the 26 irregular satel- 
lites known at the end of 2005 for Saturn, which means we 
excluded S/2004 S19 and the S/2006 satellites. Even with 
this setup, RADAU's run required about 6 months of com- 
putational time: as a comparison, HJS's run based on full 
Model 2 took approximately a month. 

The averaged (proper) orbital elements computed in Mod- 
els 1 and 2 are given in tables [l] and [2] together with the 
minimum and maximum values reached during the simula- 
tions. The same elements are visually displayed in fig.[I]and 
[2] At first sight, the mean elements obtained with Models 1 
and 2 appear to be approximately the same. However, there 
are some relevant differences between the two sets: this is 
the case of Tarvos and S/2004 S18 and of Kiviuq and Ijiraq, 
whose radial ordering results inverted. The perturbations by 
Titan and Iapetus significantly affect the secular evolution 
of the satellite system. We will show in section [3] that the 
changes are more profund that those shortly described here. 
The differences in the mean elements computed with 
RADAU and HJS on the same model (Model 2) are shown in 
the top panels of fig. [l] and [2] and are significantly less impor- 
tant: they appear to be due to the presence of chaos in the 
system. Finally, in fig. [3] we confronted the mean elements 
computed from Model 2 with the average orbital elements 
available on the JPL Solar System Dynamics websit^] In 
this case the differences between the two dataset are more 
marked, probably depending on the different integration- 
averaging time used to compute them. 

In the following discussions, when not stated differently, we 
will always be implicitly referring to the mean elements com- 
puted with Model 2 since it better represents the real dy- 
namical system. We would like to stress that the validity of 
our mean elements is proved over 10 8 years and only through 
the analysis of the satellites' secular evolution we can be able 
to assess if they can be meaningful on longer timescales. 
As a final remark, the two major gaps in the radial distri- 
bution of Saturn's irregular satellites, the first centred at 
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Figure 4. Values of the r\ parameters for the irregular satellites 
of Saturn. From top to bottom the plots show the r) values for 
eccentricity, inclination and semimajor axis respectively. The re- 
sults concern Models 1 and 2 and the additional simulation based 
on Model 1 where strict double precision (64 bit) was adopted 
in the computations instead of the standard extended precision 
(80 bit). The change in the numerical precision has major effects 
on various satellites which, as a consequence, evolve on chaotic 
orbits. 



3 Preliminary tests on a timespan of several 10 6 years, performed 
with a new library (on development) which allows the treatment 
of such effects in HJS algorithm, show no major discrepancies in 
the computed mean orbital elements and no qualitative differ- 
ences on the chaos measures we will describe in section [3] thus 
supporting the analysis and the results reported in the present 
p aper. 

4 http://ssd.jpl.nasa.gov/7sat_elein 



Phoebe and extending from 11.22 x 10 6 km to 14.96 x 10 6 
km and the second located between 20.94 x 10 6 km and 
22.44 x 10 6 km, will be discussed respectively in section [4] 
and [3] 
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Figure 5. Secular evolution of Ijiraq's longitude of pericenter in (from top left, clockwise direction) Model 2 computed with HJS 
algorithm, Model 2 computed with RADAU algorithm, Model 1, Model 1 with strict (64 bit) double precision. Angles are expressed 
in degrees and time in years. In all cases computed with HJS code, Ijiraq's longitude of pericenter is into a stable librational regime, 
indicating that the satellite is locked into a Kozai resonance with the Sun. In the case computed with RADAU, the satellite is initially 
trapped into the same resonance but it breaks the resonant regime during the last 2 X 10 7 years of the simulation. Escape from the 
resonance is not permanent, since the satellite is temporarily captured at least other two times before the end of the simulation. Ijiraq's 
behaviour is probably due to the perturbations of Titan and Iapctus. 



3 EVALUATION OF THE DYNAMICAL 
EVOLUTION 

The mean orbital elements given in the previous section give 
a global view of the dynamical evolution of the satellite sys- 
tem. However, to better understand the features of the mean 
elements and of the differences observed between the differ- 
ent models we need to have a better insight on the individual 
secular evolution of the satellites. We start our analysis by 
looking for hints of chaotic or resonant behaviour. We adopt 



a variant of the mean actions criterion (see Morbidelli ( 2002 1 
and references therein) introduced by Cuk & Burns (2004 1 
to study the dynamics of irregular satellites. This modified 
criterion is based on the computation of the r\ parameter, 
defined as 



10 

i=i 
where 



(e) ■ — — 



and 



j(T/W) 



edt 



(1) 



(2) 



(i-l)(T/10) 



<e) = 



T 



e dt. 



(3) 



The r\ parameter is a measure of the temporal uniformity 
of the secular evolution of the eccentricity. As shown in eq. 
[l] [2] and [3] its value is estimated by summing the differ- 
ences between the mean values (e) • of eccentricity computed 
on a fixed number of time intervals (10 in equations [I] and 
[2| and the mean value (e) on the whole timespan covered 
by the simulations. If the motion of the considered body 
is regular or quasi-periodic, the mean value of eccentricity 
should change little in time and the 77 parameter should ap- 
proach zero. On the contrary, if the motion is resonant or 
chaotic, the mean value of eccentricity would depend on the 
time interval over which it is computed and, therefore, r\ 
would assume increasing values. Obviously, the reliability of 
this method depends on the relationship between the num- 
ber and extension of the time intervals considered and the 
frequencies of motion: if the timescale of the variations is 
shorter than the length of each time interval, the averaging 
process can mask the irregular behaviour of the motion and 
output an r\ value lower than the one really describing the 
satellite dynamics. On the contrary, a satellite having a reg- 
ular behaviour but whose orbital elements oscillate with a 
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Figure 6. Evolution of Ijiraq's longitude of pericenter and eccen- 
tricity in Model 1 simulated with HJS algorithm using strict (64 
bit) double precision. The longitude of pericenter and the eccen- 
tricity have limited non-periodic changes which may be ascribed 
to the perturbations of Jupiter since Titan and Iapetus are not 
included in the model. 



period longer than the considered time intervals may spuri- 
ously appear chaotic. There is no a priori prescription con- 
cerning the number of time intervals to be used: after a set 
of numerical experiments, we opted for the value (10 time 
intervals) suggested in the original work by Cuk & Burns 
(20041. 



In the following we extend the idea of Cuk & Burns ( 2004 1 



by applying a similar analysis also to inclination and semi- 
major axis, thus introducing r\ a and rji as: 



10 
3 = 1 

and 



3=1 



(H-W)' 

(a) 2 
« 2 



(4) 



(5) 



with corresponding definitions of the quantities involved. 
We present the r\ values in fig. |4j where we compare the 
output of Model 1, of Model 2 and of an alternative version 
of Model 1 where a different precision in the floating point 
calculations have been used. We we will shortly explain the 
reason of this duplication. 

The coefficient showing the largest variations is rj s (see fig. 
[4] top panel) and the extent of the dynamical variations for 
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Figure 7. Evolution of Ijiraq's longitude of pericenter, eccen- 
tricity and inclination (from top to bottom) in Model 2 inte- 
grated with RADAU. The dynamical evolution of the satellite 
has a sudden change around 6 X 10 7 years leading to a progres- 
sive increase in the eccentricity oscillations. Also the libration 
amplitude slightly grows until the longitude of pericenter starts 
to circulate at 8 X 10 7 years. Since fig. [H] proved that the effects 
of the giant planets are quite limited after the onset of the Kozai 
resonance with the Sun, Titan and Iapetus are responsible for the 
perturbed evolution of the satellite. 



each satellite can be roughly evaluated by the square root 
of ry e . These variations are on average around 5 — 10% with 
peak values of about 30%. The corresponding variations in 
semimajor axis and inclination, as derived from the rj a and 
rji, are on average more limited ranging from 1% to 2%. Since 
the values of rj e are not significantly correlated to those of 
7] a and r]i, our guess is that the different coefficients are in- 
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Figure 8. Secular evolution of Kiviuq's longitude of pericenter. From top left in clockwise direction we plot the outcome of Model 2 
computed with HJS algorithm, Model 2 computed with RADAU algorithm, Model 1, Model 1 with strict (64 bit) double precision. 
Angles are expressed in degrees and time in years. Transitions between circulation and libration appear in both the dynamical models 
and independently from the numerical codes. Librations in Model 2 - HJS occur between 6 X 10 7 and 7 X 10 7 years and are extremely 
short-lived (the angle mainly circulates). The librations are also evenly distributed around 90° and 270°. The same model computed 
with RADAU shows long-lived librational periods centred at 270°. Also in Model 1 the librations are long-lived but they concentrated 
around 90° . In Model 1 (bottom right panel) at the end of the simulation the satellite enters a libration phase which lasts about 2 X 10 7 
years. The differences between the four cases argue for the chaotic nature of Kiviuq's dynamical evolution. Titan and Iapetus intervene 
altering the dynamical history of the satellite. 



dicators of different dynamical effects. 

We used a Model 1 characterised by a different numerical 
precision in the calculations to unmask the presence of chaos 
in the system. The standard set-up on x86 machines for dou- 
ble precision computing is the so called extended precision, 
where the floating point numbers are stored in double pre- 
cision variables (64 bit) but a higher precision (80 bit) is 
employed during the computations. This is also the set-up 
we originally used in Model 1, which in the graphs is labelled 
"Model 1 (80 bit)". We ran an additional simulation forcing 
the computer to strictly use double precision, thus employ- 
ing 64 bits also during the computations. The results of this 
simulation are labelled as "Model 1 (64 bit)" . Hereinafter, 
unless differently stated, when we refer generically to Model 
1 we intend the 80 bit case. From a numerical point of view, 
the change in computing precision should in principle affect 
only the last few digits of each floating point number, since 
the values are stored in 64 bit variables in both cases. Reg- 
ular motions should be slightly affected by the change while 
chaotic motions should show divergent trajectories. By com- 
paring the different values of the rj parameters for "Model 
1 (64 bit)" and "Model 1 (80 bit)" we find some satellites 



with drastically different values, a clear indication of chaotic 
motion. 

After this preliminary global study of the satellite system, 
we have performed a more detailed analysis of selected ob- 
jects whose behaviours have been already investigated in 
previous papers. Carruba et al. (20021, Nesvorny et al. 



20031) and Cuk & Burns ( 2004 1 reported four cases of res- 



onant motion between the irregular satellites of Saturn: 

• Kiviuq and Ijiraq being in the Kozai regime with the 
Sun; 

• Siarnaq and Paaliaq, whose pericenters appeared tidally 
locked to the Sun. 

The authors computed the orbital evolution of these satel- 
lites with a modified version of Swift N-Body code, where 
the planets were integrated in the heliocentric reference 
frame while the irregular satellites in planetocentric frames. 
This dynamical structure is similar to that of our Model 1. 
We extracted from our simulations the data concerning the 



same objects studied in Carruba et al. (20021, Nesvorny et 



al] ( |2003| ) and |Cuk fc Burns | ( |2004[ )~and we analysed their 
behaviour. The motion of Ijiraq as from our Models 1 and 2 
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Figure 9. Comparison between the evolution of semimajor axis, eccentricity and inclination (from top to bottom) of Paaliaq (left 
column) and Siarnaq (right column) in Model 2 integrated with HJS. Both satellites suffer a slight inward displacement, limited to about 
~ 1 — 2%, yet Paaliaq's eccentricity and inclination show a regular behaviour wher eas Siar naq 's ones have a coupled variation at about 
1 X 10 7 years and 5 X 10 7 years. This behaviour is reproduced in both Models (see fig, |lO| and |lf | for further details on Siarnaq's evolution). 



integrated with HJS (see fig. [5] first three panels from top 
left in counterclockwise direction) appears to evolve in a sta- 
ble Kozai regime with the Sun. The longitude of pericenter 
librates around 90° with an amplitude of ±30° . The simula- 
tions based on Model 2 and Model 1 (80 bit precision) show 
a regular secular behaviour of the orbital elements while the 
simulation based on Model 1 (64 bit precision) shows peri- 
ods in which the range of variation of the eccentricity shrinks 
by a few percent in correspondence to a similar behaviour of 
the longitude of pericenter (see fig.[6]for details). The output 
obtained from Model 2 with RADAU algorithm (fig.[7| is in- 
stead significantly different. Ijiraq is in a stable Kozai regime 



for the first 6 x 10 7 years, then it experiences a change in the 
secular behaviour of both eccentricity and inclination and, 
finally, the longitude of pericenter va circulates for about 
1.3 x 10 7 years. From then on, zu alternates phases of circu- 
lation and libration around both 90° and 270°. 
Kiviuq shows an even more complex behaviour (see fig. Isb. 
All the simulations (Models 1 and 2 and both HJS and 
RADAU algorithms) predict transitions between periods 
of circulation and libration of zu. The centre of libration 
changes depending on the numerical algorithm: the secular 
evolution computed with RADAU algorithm show a preva- 
lence of librational phases around 270° . In the case of HJS 
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Figure 10. Secular evolution of the eccentricity of Siarnaq. From top left, clockwise direction, the plots show the outcome of Model 2 
- HJS algorithm, Model 2 - RADAU algorithm, Model 1 - HJS with standard precision and Model 1 - HJS with strict (64 bits) double 
precision. Siarnaq's eccentricity clearly shows evidence of chaotic evolution in Model 2 and Model 1 (64 bit). In Model 1 (80 bit) Siarnaq's 
eccentricity appears quasi-stationary at all timescales. 



algorithm the dominant libration mode is around 90°. While 
the libration cycles were all characterised by an amplitude 
of about ±30°, their durations vary from case to case. In 
Model 2 computed with HJS w mostly circulates with only 
a few short-lived librational periods. In Model 1 (80 bit) 
the satellite enters the Kozai regime with the Sun only af- 
ter 8 x 10 7 years, while Model 1 (64 bit) shows a behaviour 
more similar to that of Model 2, even if with longer-lived 
librational phases. Model 2 computed with RADAU shows a 
behaviour similar to that of Model 1. As for Ijiraq, such dif- 
ferent dynamical evolutions confirm the chaotic behaviour 
of the satellite orbit. 

Ijiraq and Kiviuq are examples of a different balancing be- 
tween the perturbations acting on the trajectories of these 
satellites. Ijiraq appears dominated by the Kozai resonance 
caused by the Sun's influence while Titan, Iapetus and the 
other planets play only a minor role. On the opposite side, 
Kiviuq is significantly perturbed by Titan, Iapetus and the 
other planets which prevent its settling into a stable Kozai 
regime with the Sun and cause its chaotic evolution. 
For Siarnaq and Paaliaq our updated dynamical model do 
not confirm the resonant motions found in previous publica- 
tions. In presence of Titan and Iapetus (Model 2) both the 
satellites show evidence of a limited but systematic varia- 
tion of their semimajor axes. They monotonically migrate 
inwards by ~ 2% during the timespan covered by our sim- 



ulations (see fig. |9j. Paaliaq has a regular evolution of the 
eccentricity and inclination (fig. [9] left column), Siarnaq (fig. 
[9] right column) during its radial migration presents phases 
of chaotic variations of both eccentricity and inclination (see 
fig. |10| and |ll| for details). We cannot rule out major changes 
in semimajor axis (possibly periodic?) on longer timescales. 
Paaliaq, during its orbital evolution, might enter more dy- 
namically perturbed regions (see fig. [26] for further details). 
The comparison of our results with previously published pa- 
pers suggest that the three-body approximation adopted in 
previous analytical works was not accurate enough to be 
used as a reference model for the dynamical evolution of 
Saturn's irregular satellites. An additional feature arguing 
against a simplified three-body approximation is the fol- 
lowing one. When we illustrate the secular evolution of the 
satellites in the e — i plane, in several cases we obtained a 
thick arc (see fig. [12] for details) having opposite orientation 
for prograde and retrograde satellites. It indicates a global 



anticorrelation of the two orbital elements (see fig. 13 for 



details) and it appears more frequently among the satellites 
integrated with Model 2. This anticorrelation of eccentric- 
ity and inclination is a characteristic feature of the Kozai 
regime. The arc-like feature is in fact manifest in the two 
previously discussed cases: Kiviuq and Ijiraq. By inspecting 
the dynamical histories of all the other satellites, we found 
the following with the same feature: Paaliaq, Skathi, Al- 
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Figure 11. Secular evolution of the inclination of Siarnaq. From top left, clockwise direction, the plots show the outcome of Model 
2 - HJS algorithm, Model 2 - RADAU algorithm, Model 1 - HJS with standard precision and Model 1 - HJS with strict (64 bits) 
double precision. While Paaliaq radial displacement is coupled to a regular evolution of the other orbital elements (see fig. [9j , Siarnaq's 
eccentricity (see fig. |10| l and inclination show clearly the presence of chaotic features in Model 2 and Model 1 (64 bit). In Model 1 (80 
bit) Siarnaq's inclination appears quasi-stationary at all timescales. 



biorix, Erriapo, Siarnaq, Tarvos, Narvi, Bebhionn (S/2004 
Sll), Bestla (S/2004 S18), Hyrokkin (S/2004 S19). Three 
of these satellites have an inclination close to the critical 
one leading to Kozai cycles. Farbauti (S/2004 S9), Kari 
(S/2006 S2), S/2006 S3 and Surtur (S/2006 S7) show a more 
dispersed arc-like feature suggesting that external gravita- 
tional perturbations by the massive satellites and planets 
interfere with the Kozai regime due to solar gravitational 
effects. This hypothesis is confirmed by the fact that in a 
few cases the arc-like feature is present only in Model 2 
were Titan and Iapetus are included. 

Additional differences in the chaotic evolution of some satel- 
lites can be ascribed to the presence of Titan and Iapetus. 
According to our simulations 23 satellites show a chaotic 
behaviour. Aside from the previously mentioned Kiviuq, Iji- 
raq, Paaliaq and Siarnaq, the list is completed by: Albiorix, 
Erriapo, Tarvos, Narvi, Thrymr, Ymir, Mundilfari, Skathi, 
S/2004 S10, Bebhionn (S/2004 Sll), S/2004 S12, S/2004 
S13, S/2004 S18, S/2006 S2, S/2006 S3, S/2004 S4, S/2006 
S5, S/2006 S7 and S/2006 S8 (see fig.[I3|. Some of them are 
chaotic only in one of the two dynamical Models (e.g. fig. 



15 fig. 16 and fig. 17 1 while others in both the Models (see 
fig. 18 as an example). The most appealing interpretation is 



that Titan and Iapetus perturb the dynamical evolution of 
the irregular satellites either by stabilising otherwise chaotic 



feature is that in most cases the chaotic features appeared 
in the secular evolution of the semimajor axis with eccen- 
tricity and inclination exhibiting a regular, quasi-periodic 
evolution. As a consequence, the present structure of the 
outer Saturnian system could be not representative of the 
primordial one. We will further explore this issue in section 
® 

The presence of chaos in the dynamical evolution of the 
irregular satellites appears to be the major driver of the dif- 
ferences between the outcome of Model 1 (80 and 64 bit) and 
2 (HJS and MERCURY RADAU algorithms). The chaotic 
nature of the trajectories is probably also at the origin of 
the differences (e.g. the alternation between resonant and 
not resonant phases in Kiviuq's evolution) we noticed be- 
tween our simulations based on Model 1 and those published 



in the literature by Carruba et al. (20021, Nesvorny et al 



( 2003 1 and Cuk & Burns ( 2004 1, which were based on a sim 



orbits (fig. 16 1 or causing chaos (fig. 17 1. Another interesting 



ilar dynamical scheme. The interplay between the inclusion 
of Titan's and Iapetus' gravitational perturbations and the 
presence of chaos can finally be invoked to explain the dif- 
ferences between the results obtained with Models 1 and 2. 
To explore in more detail the effects of Jupiter, Titan and 
Iapetus on the dynamics of irregular satellites we performed 
two additional sets of simulations where we sampled with 
a large number of test particles the phase space populated 
by the satellites. We distributed 100 test particles in be- 
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Figure 12. Anticorrelation of inclination and eccentricity (see also fig. |13| in the dynamical evolution of (clockwise from top left) Ijiraq, 
Paaliaq, Bestla (S/2004 S18) and Narvi. These features, observed in both Models 1 and 2, might imply that these satellites were prevented 
from entering the Kozai cycle with the Sun by the combined gravitational perturbations of the other giant planets. It is noteworthy that 
the only satellites in a full Kozai regime are the two innermost ones, for which the external perturbations are less effective. 



tween 10.47 x 10 6 km and 25.43 x 10 6 km on initially cir- 
cular orbits and integrated their trajectories for 10 6 years. 
10 simulations have been performed each with a different 
value of the initial inclination of the particles. We con- 
sidered i = 0°, 15°, 30°, 45°, 60° for prograde orbits, and 
i = 120°, 135°, 150°, 165°, 180° for retrograde orbits. To 
compute the orbital evolution we used the HJS N-Body 
code and the dynamical structure of Models 1 and 2. The 
initial conditions for the massive perturbing bodies (Sun, 
giant planets and major satellites) and the timesteps were 
the same used in our previous simulations. 
To analyse the output data we looked at the mean elements 
and the r\ parameters. In fig. |19| and|20|we show some exam- 
ples of the outcome where the influence of Titan and Iape- 
tus is manifest while comparing the outcome of Model 1 vs. 
Model 2. For particles in the inclination range [30° — 60°] 
and [120° — 150°] the values of the r/e are significantly dif- 
ferent. In addition, we verified the leading role of Jupiter in 
perturbing the system. We computed the test particle or- 
bits switching off the gravitational attraction of the planet. 
The dynamically perturbed regions shown fig. [2l] disappear 
or are reduced (Titan and Iapetus are still effective) when 
Jupiter's pull is switched off. 

By inspecting the mean orbital elements of irregular satel- 
lites in fig. [2l] one may notice as those orbiting close to 
the perturbed regions have a "stratification" in eccentricity 



typical of a resonant population. It is possible that the dy- 
namical features we observe today are not primordial but 
a consequence of a secular interaction with Jupiter. One of 
these features is the gap between about 20.94 x 10° km and 
22.44 x 10 6 km in the radial distribution of the mean ele- 
ments of the irregular satellites. This gap was particularly 
evident till June 2006, since no known satellite populated 
that region. With the recent announcement of 9 additional 
members of Saturn's irregular satellites, two objects are now 
known to cross this region: S/2006 S2 and S/2006 S3. Both 
these bodies have their mean semimajor axis falling near to 
perturbed regions. According to our simulations, the semi- 
major axis of both satellites show irregular jumps (see fig. 
|22[ ) typical of chaotic behaviours. S/2006 S2 has phases of 
regular evolution alternated with periods of cha os lasting 
from a few 10 4 up to about 10 s years (see fig. 



22 



left plot). 

In the case of S/2006 S3, we observed a few large jumps 
in semimajor axis over the timespan covered by our simu- 



lations (see fig. 22 right plot). The evolutions of both the 



eccentricity and the inclination of the satellites appear more 
regular, with long period modulations and beats. 
We applied the frequency analysis in order to identify the 
source of the perturbations leading to the chaotic evolution 
of the two satellites. By analysing the h and k non singu- 
lar variables defined with respect to the planet, we found 
among the various frequencies two with period of 900 years 
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Figure 13. Phase displacement of inclination and eccentricity in the dynamical evolution of (clockwise from top left) Ijiraq, Paaliaq, 
Bestla (S/2004 S18) and Narvi. Both eccentricity and inclination have been rescaled and shifted to show variations of the same order of 
magnitude. The plots referring to prograde satellites show the relation between eccentricity and inclination, those referring to retrograde 
satellites show the relation between eccentricity and the inclination supplemental angle (i.e. 180° — i). Bestla's evolution show alternating 
phases of correlation and anticorrelation between eccentricity and inclination: anticorrelation phases generally last twice than correlation 
phases. This effect is at the basis of Bestla's thicker arc in fig. |12| The evolutions of Bestla's inclination and eccentricity are strongly 
coupled (sec fig. |27] for details) and this coupling interferes with solar perturbations causing Kozai regime. 
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Figure 14. Chaotic and resonant irregular satellites of Saturn's system. The square symbols identifies the irregular satellites showing 
resonant or chaotic features in their dynamical evolution. The data are presented in the a — e (left plot) and a — i (right plot) planes. 



and 1800 year respectively. These values are close to the 
Great Inequality period (~ 883 years) of the almost res- 
onance betwen Jupiter and Saturn. This is possibly one 
source of the chaotic behaviour of the satellites. In addi- 
tion, the inspection of the upper plot in fig. [20] shows that 



the rj a parameters of the test particles populating this ra- 
dial region increase when Titan and Iapetus are included. 
By comparing the frequencies of motion of the two irregular 
satellites with those of Titan and Iapetus we find an addi- 
tional commensurability. The frequencies 7.1 x 10 -4 yr~ x 
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Figure 15. Secular evolution of the inclination of S/2006 S7 in Model 2 (left panel) and Model 1 computed with strict (64 bits) double 
precision (right panel). While the first case show a regular behaviour, the last one presents an evident chaotic jump in inclination. This 
jump may be related to the entry of the satellite into a chaotic region. This event does not occur in the Model 1 computed with standard 
precision, probably due to the different numerical setup. 
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Figure 16. Secular evolution of the semimajor axis of Mundilfari. From top left, clockwise direction, the plots show the outcome of 
Model 2 - HJS algorithm, Model 2 - RADAU algorithm, Model f - HJS with standard precision and Model 1 - HJS with strict (64 bits) 
double precision. The simulations based on Model 2 show a regular behaviour on the whole timespan, independently of the employed 
algorithm. Both cases based on Model 1 show instead an irregular evolution diagnostic of chaotic behaviour. Titan and Iapetus act like 
orbital stabilisers for Mundilfari. 



and 7.4 x 10~ 4 yr' 1 that are present in the power spec- 
trum of S/2006 S2 and S/2006 S6, respectively, are about 
twice the frequency 3.4 x 10 -4 yr~ x in the power spectrum 
of Titan. The cumulative effects of the Great Inequality and 
of Titan and Iapetus lead to the large values of the r\ pa- 



rameters in fig. [4] and |20| The chaotic evolution of the two 
satellites does not lead to destabilisation in the timespan of 
our integration (10 8 years), however longer simulations are 
needed to test the long term stability. It is possible that the 
irregular behaviour takes the two satellites into other chaotic 
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Figure 17. Secular evolution of the semimajor axis of Erriapo. From top left, clockwise direction, the plots show the outcome of Model 
2 - HJS algorithm, Model 2 - RADAU algorithm, Model 1 - HJS with standard precision and Model 1 - HJS with strict (64 bits) double 
precision. Contrary to the case showed in fig. |16| Titan and Iapetus cause an irregular evolution of Erriapo. In Model 1 the satellite has 
a regular behaviour. 
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Figure 18. Secular evolution of the eccentricity of S/2006 S5 in Model 2 computed with HJS (left panel) and Model 1 computed using 
standard precision (right panel). The chaotic nature of the satellite orbital motion can be easily inferred by the secular behaviour of its 
eccentricity. Similar conclusions can be drawn from the inspection of semimajor axis and inclination. 



regions ultimately leading to their expulsion from the sys- 
tem. It is also possible that, during their chaotic wandering, 
they cross the paths of more massive satellites and be colli- 
sionally removed. This could have been the fate of possible 
other satellites which originally populated the region encom- 
passed between 20.94 x 10 6 km and 22.44 x 10 6 km. 



4 EVALUATION OF THE COLLISIONAL 
EVOLUTION 

To understand the present orbital structure of Saturn's 
satellite system and how it evolved from the primordial one 
we have to investigate the collisional evolution within the 
system. Impacts between satellites, in fact, may have re- 
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Figure 19. Evolution of the r\ parameters with semimajor axis for prograde test particles in Model 1 (left column) and Model 2 (right 
column). From top to bottom, the plots show the values of r? e , r\ a respectively. Values of r\ near indicate quasi periodic motion, 
higher values indicate increasingly chaotic or resonant behaviour. Distances are expressed in 10 6 km. 



moved smaller bodies and fragmented the larger ones. Mi- 
nor bodies in heliocentric orbits like comets and Centaurs 
may have also contributed to the system shaping by colliding 
with the satellites as addressed by |Nesvorny et ah] |2004). 
At present, however, such events are not frequent because 
of the reduced flux of minor bodies and the small sizes of 
irregular satellites (Zahnle et al. 2003, Ncsvo rny et al.|2004[ ). 
The last detailed evaluation of the collisional probabilities 
between the irregular satellites around the giant planets was 
the one performed by |Nesvorny et al.| (j2003 ) , which showed 
that the probabilities of collisions between pairs of satellites 
were rather low and practically unimportant. The computed 
average collisional lifetimes were tens to hundreds of times 



longer than Solar System's lifetime. The only notable ex- 
ceptions were those pairs involving one of the big irregu- 
lar satellites (e.g. Himalia, Phoebe, etc.) in the systems. In 
the Saturn system, Phoebe is between one and two order 
of magnitudes more active than any other satellite. The au- 
thors computed a cumulat ive number of collisio ns between 
6 and 7 in 4.5 x 10 9 years ( |Nesvorny et al.|2003 l. 
However, at the time of the publication of the work by 



Nesvorny et al. ( 2003 1 only 13 of the 35 currently known ir- 



regular satellites of Saturn had been discovered. We extend 
here their analysis taking advantage of the larger number 
of known bodies and of the improved orbital data. Using 
the mean orbital elements computed with Model 2 we have 
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Figure 20. Evolution of the r\ parameters with semimajor axis for retrograde test particles in Model 1 (left column) and Model 2 (right 
column). From top to bottom, the plots show the values of r? e , r\ a respectively. Values of r\ near indicate quasi periodic motion, 
higher values indicate increasingly chaotic or resonant behaviour. Distances are expressed in 10 6 km. 



estimated the collisional probabilities using the approach de- 



scribed by Kcssler ( 1981 1. Since in the scenario described by 



the Nice Model the Late Heavy Bombardment (LHB) rep- 
resents a lower limit for the capture epoch of the irregular 



satellites ( Gomes et al. 2005 Tsiganis et al. 2005 ) , we consid- 



ered a time interval for the collisional evolution of 3.5 x 10 9 
years, making the conservative assumption that the LHB 
took place after about 10 9 years since the beginning of Solar 
System formation. Since the collisional probability depends 
linearly on time, our results can be immediately extended 
to longer timescales. 

The results of our computations (see fig. 



23 ) confirm 



that the only pairs of satellites with high probability 



of collisions involve Phoebe. The satellite pairs with the 
higher number of collisions over the considered timespan 
are Phoebe-Kiviuq (0.7126), Phoebe-Ijiraq (0.7099) and 
Phoebe-Thrymr (0.6436). The remaining pairs involving 
Phoebe have a number of collisions between 0.1 and 0.35 



(see fig. 23 line/column 3), with the highest values associ- 
ated with the prograde satellites Paaliaq, Siarnaq, Tarvos, 
Albiorix, Erriapo and Bebhionn (S/2004 Sll). All the other 
satellite pairs, due also to their small radii, have negligi- 
ble (< 10~ 2 ) collisional probabilities. The predicted total 
number of collisions in Saturn's irregular satellite system, 
obtained by summing over all the possible pairs, is of about 
12 collisions over the considered timespan. Half of these im- 
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Figure 21. Location of the irregular satellites in the a — e (upper plot) and a — i (lower plot) planes respect to the unstable regions 
identified by the simulations computing the evolution of regularly sampled test particles. The unstable regions, perturbed by Jupiter, are 
marked by dashed and dotted lines. Dotted lines indicate regions which are excited for all values of inclination. Dashed lines are related 
to those regions which get perturbed for high absolute inclination values (i.e. i ~ 150°). 
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Figure 22. Secular evolution of the semimajor axis of S/2006 S2 Kari (detail, left graph) and S/2006 S3 (detail, right graph) in model 
2. While on different timescales and with different intensities, both the satellites have a chaotic evolution possibly due to Jupiter's 
perturbations through the Great Inequality with Saturn. 



pacts involve Phoebe. This is probably at the origin of the 
gap centred at Phoebe and radially extending from about 
11.22 x 10 6 km to about 14.96 x 10 6 km from Saturn (i.e. be- 
tween Ijiraq's and Paaliaq's mean orbits) for both prograde 
and retrograde satellites. To further confirm this hypothesis 
we evaluated the impact probability for a cloud of test par- 
ticles populating this region. 

We filled with test bodies a volume in the phase space de- 
fined in the following way: 

• 10.47 x 10 6 km < a sC 14.96 x 10 6 km 

• < e sC 0.9 

• 0° < i $J 50° for prograde orbits 

• 130° ^ i < 180° for retrograde orbits 

The sampling stepsize were 5a = 4.488 x 10 4 km, 8e — 0.1 
and Si — 2°, for a total of 2860 prograde orbits and 2860 
retrograde ones. A radius of 3 km has been adopted for the 
test particles to compute the cross sections. The results are 
presented in the colour maps of fig.[24]for the prograde cases 
and [25] for the retrograde ones. Our results show that low- 
est collision probabilities (of the order of 0.1 collisions) are 



related to orbits with high values of eccentricity and inclina- 
tion (e > 0.5 and i > 25° or i < 155°) as illustrated in the 
top right quadrants of each plot of fig. [24] and the bottom 
right quadrants of each plot of fig. |25| Prograde low inclina- 
tion orbits with 0° < i ^ 10° and e < 0.5 have more than 
4 collisions over the given timespan (see the bottom left 
quadrants of each plot of fig. 241. Retrograde orbits with 
174° < i < 180° have at least one collision for any value of 



eccentricity (see the top part of each plot of fig. 25 1. The 
number of collisions with Phoebe is generally higher for the 
prograde test particles (1 — 5 collisions) compared to the 
retrograde orbits (0.5 — 1). These results support the ex- 
istence of a Phoebe's gap caused by collisional erosion. It 
cannot be due to dynamical clearing mechanisms since, ac- 
cording to an additional set of simulations, we show that the 
region around the Phoebe's gap is stable. 100 test particles 
are uniformly distributed across the Phoebe's gap with the 
following criteria: 

• semimajor axis from 10.47 x 10 6 km to 14.96 x 10 6 km 

• radial spacing between the particles of 4.488 x 10 4 km 

• initially circular orbits 
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against the hypothesis of a Phoebe family. The existence 
of Phoebe's family has been a controversial subject since 



its proposition by Gladman et al. (2001 1. Its existence was 
guessed on the close values of inclination of Phoebe and 
other retrograde satellites. The dynamical inconsistency of 



this criterion has been pointed out by Nesvorny et al. ( 2003 1, 



who showed that the velocity dispersion required to relate 
the retrograde satellites to Phoebe would be too high to be 
accounted as realistic in the context of the actual knowl- 
edge of fragmentation and disruption processes. Our results 
suggest also that if a breakup event involved Phoebe, the 
fragments would have been ejected within the Phoebe's gap 
for realistic ejection velocities. As a consequence, they would 
have been removed by impacting on Phoebe. 



Figure 23. Colour map of the collisional probabilities between 
the 35 presently known irregular satellites of Saturn. For each 
pair of satellites the probability of collision is evaluated in terms 
of the number of impacts during the system lifetime (3.5 X 10 9 
years). The colour coding used to illustrate the impact probabil- 
ities ranges from yellow (highest values) to blue (lowest values). 
Black is used when the orbits do not intersect. The only pairs 
of satellites for which the impact probability is not negligible al- 
ways include Phoebe as it can be argued by the fact that only 
column/line 3 (Phoebe) is characterised by bright colours. 



• initial inclinations set to fixed values equal 
to 0°, 15°, 30°, 45°, 60° for prograde orbits and 
120°, 135°, 150°, 165°, 180° for retrograde orbits. 

The orbits have been integrated for 10 6 years and for each 



5 EXISTENCE OF COLLISIONAL FAMILIES 

The identification of possible dynamical families between the 
irregular satellites of the giant planets had been a common 
task to all the studies performed on the subject. The exis- 
tence of collisional families could in principle be explained by 
invoking the effects of impacts between pairs of satellites and 
between satellites and bodies on heliocentric orbits. The im- 
pact rate between satellite pairs is low even on timescales of 
the order of the Solar System age, with the only exception of 
impacts amongst the most massive irregular satellites. The 
gravitational interactions between the giant planets and the 
planetesimals in the early Solar System may have pushed 
some of them in planet-crossing orbits. This process is still 
active at the present time and Centaurs may cross the Hill's 



of them a value of r\ has been computed (see Fig. 26 1. The r\ sphere of the planets. However, Zahnle et al. (20031 showed 



values indicate that the perturbations of Titan and Iapetus 
are negligible and that test particles positioned in regions of 
high collision probability (i.e. for i 10° and i ^ 174°) are 
not dynamically unstable. 

The collisional origin of the Phoebe's gap is also confirmed 
by observational data showing that the irregular satellites 
moving closer to Phoebe are those located in regions of the 
phase space where the impact probability with Phoebe is 
lower. In addition, the images of Phoebe taken by ISS on- 
board the Cassini spacecraft revealed a strongly cratered 
surface, with a continuous crater size distribution ranging 
from about 50 m, a lower value imposed by the resolution of 
ISS images, to an upper limit of about 100 km comparable 
to the dimension of the satellite. A highly cratered surface 



that the present flux of bodies, combined with the small size 
of irregular satellites, is unable to supply an adequate im- 
pact rate. If collisions between planetesimals and satellites 
are responsible for the formation of families, these events 
should date back sometime between the formation of the gi- 
ant planets and the Late Heavy Bombardment. 
The possible existence of dynamical families in the Saturn 
satellite system has been explored by using different ap- 
proaches. Photometric comparisons has been exploited by 



Grav et al. (20031; Grav & Holman (20041; Buratti et al 



( |2005[ ) and were limited to a few bright objects. Dynamical 



methods have been used by Grav et al. (20031; Nesvorny 



et al. (20031; Grav & Holman (20041 but on the limited 



was predicted by Nesvorny et al. ( 2003 1 , who also suggested 



that the vast majority of Phoebe's craters should be due to 
either impacts with other irregular satellites (mainly pro- 
grade ones) or be the result of a past intense flux of bodies 
crossing Saturn's orbit. Comets give a negligible contribu- 
tion having a frequency o f collision with Pho ebe of about 
1 impact every 10 9 years (Zahnle et al. 20031. Our results 



sample (about one third of the presently known population) 
of irregular satellites available at that time. These methods 
aimed to identify those satellites which could have origi- 
nated from a common parent body following one or more 
breakup events. The identification was based on the evalua- 
tion through Gauss equations of the dispersion in orbital el- 
ement space due to the collisional ejection velocities. In this 
paper we apply the Hierarchical Clustering Method (here- 



confirm those of Nesvorny et al. (20031 showing that indeed inafter HCM) described in Zappala et al. (1990 19941 to 



Phoebe had a major role in shaping the structure of Sat- 
urn's irregular satellites. The existence of a primordial now 
extinct population of small irregular satellites or collisional 
shards between 11.22 x 10 6 km and 14.96 x 10 6 km from Sat- 
urn could explain in a natural way the abundance of craters 
on Phoebe's surface. 

Phoebe's sweeping effect appears to have another major con- 
sequence, related to the existence of Phoebe's gap: it argues 



the irregular satellites of Saturn. HCM is a cluster-detection 
algorithm which looks for groupings within a population of 
minor bodies with small nearest-neighbour distances in or- 
bital element space. These distances are translated into dif- 
ferences in orbital velocities via Gauss equations and the 
membership to a cluster or family is defined by giving a 
limiting velocity difference (cutoff ). |Nesvorny et al.| ( |2003[ ) 
adopted a cutoff velocity value of 100 m/s according to hy- 
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Figure 24. Colour map of the collisional probabilities between Phoebe and a cloud of prograde, 3 km radius test particles filling the 
region of Phoebe's gap at six values of semimajor axis: 10.47xl0 6 , 11.37xl0 6 , 12. 27 X 10 s , 13.16 X 10 6 , 14.06 xlO 6 , 14.96 X 10 6 km. The 
probabilities of collision are evaluated in terms of the number of impacts during the system lifetime (3.5 X 10 9 years) and are represented 
through a colour code going from yellow (highest values) to blue (lowest values) with black representing the case in which the orbits 
do not intersect. Note that, due to the high values of collisional probabilities reached in the prograde cases (up to 30 collisions in the 
timespan considered) we had to limit the colour range to a maximum of 10 collisions to maintain the scale readability. 



drocode models (Benz & Aspaugh||1999 1. Here we prefer to 



relax this value to 200 m/s considering the possible range of 
variability of the mean orbital elements because of dynami- 
cal effects. The results we obtained are summarised in table 
[3] and interpreted as in the following. 

By inspecting our data, we conclude that, as already ar- 



gued by Nesvorny et al. (20031, the velocity dispersion of 



prograde and retrograde satellites (about 500 m/s for pro- 
grades and more than 600 m/s for retrogrades) makes ex- 
tremely implausible that each of the two groups originated 



by a single parent body. In addition, the classification in 
dynamical groups based on the values of the orbital incli- 



nation originally proposed by Gladman et al. (20011 and 



reported by other authors (see Sheppard ( 2006[) and refer- 



ences within) is probably misleading. Prograde satellites like 
Kiviuq, Ijiraq, Siarnaq and Paaliaq do share the same incli- 
nation but, as a group, they have a velocity dispersion of 
over 450 m/s, hardly deriving from the breakup of a single 
parent body. The enlarged population of retrograde satel- 
lites we have analysed show that the clustering around a 
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Figure 25. Colour map of the collisional probabilities between Phoebe and a cloud of retrograde, 3 km radius test particles filling the 
region of Phoebe's gap at six values of semimajor axis: 10.47xl0 6 , 11.37xl0 6 , 12.27xl0 6 , 13.16 X 10 6 , 14.06 xlO 6 , 14.96 X 10 6 km. The 
probabilities of collision are evaluated in terms of the number of impacts during the system lifetime (3.5 X 10 9 years) and are represented 
through a colour code going from yellow (highest values) to blue (lowest values) with black representing the case in which the orbits do 
not intersect. 



single inclination (Gladman et al. 20011 and their associa- 
tion to Phoebe ( Gladma n et al.|2001 1 is not an indication of 
a common origin. The required velocity dispersion is in fact 
about 650 m/s. 

We found two potential dynamical families between the pro- 
grade satellites: the couple Kiviuq-Ijiraq and what we term 
as Albiorix family, composed of Albiorix, Erriapo, Tarvos 
and S/2004 Sll. The analysis of the retrograde satellites is 
more complex. There are three possible groups each com- 
posed of three satellites and two others by two satellites, 
all characterised by acceptable values of the velocity dis- 



persion (between 100 and 170 m/s). A sixth possible group 
satisfying our acceptance criterion is composed of Narvi and 
S/2004 S18, but its interpretation is quite tricky. This group 
shows a high velocity dispersion, at the upper limit of our 
range, but the dynamical evolution of both satellites is un- 
certain on a timescale of 10 9 years. The orbits of both bodies 
have the most extreme values of inclination among all retro- 
grade satellites. Our numerical experiments with test parti- 
cles showed that for such bodies the eccentricity is strongly 



coupled to the inclination (see fig. 271. In our simulations, 



initially circular orbits became highly eccentric in less than 
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Figure 26. Comparison between the r\ e parameters of prograde (top panels) and retrograde (bottom panels) test particles moving close 
to the Phoebe's gap (i.e., the region extending from 10.47 X 10 6 km to 14.96 X 10 6 km from Saturn) in Models 1 (left column) and 2 
(right column). By comparing the plots pairwise, we note that the differences increase with larger inclinations of the test particles. Those 
orbits with inclination in the range 30° — 60° for the prograde case and 120° — 150° for the retrograde one are the most influenced by the 
presence of Titan and Iapetus. This hold particularly true in the inner part of this region (i.e. up to Phoebe's orbital distance). Prograde 
test particles are less affected by the presence of Titan and Iapetus than retrograde ones. In Model 1, the retrograde test particles have 
smaller values of r] e than the prograde ones. In Model 2, the retrograde test particles show a marked increase in their rj e values for 
high inclination orbits (i.e. ~ 120°). Such orbits should likely be subject to a strongly chaotic evolution and could be short-lived on a 
timescale longer than the one (10 6 years) we considered. A noteworthy feature present in all plots is the peak at ~ 13.015 X 10 e km (the 
orbital distance of Phoebe) for both prograde and retrograde test particles. A second peak, peculiar to retrograde orbits, is located at 
~ 14.21 X 10 6 km. These peaks appear only for inclination of 30° (150° for retrograde orbits). These values bring the orbits of the test 
particles nearer to the equatorial plane of Saturn (once corrected for the giant planet's axial tilt) and thus nearer to the orbital plane of 
Titan and, to a minor extent, Iapetus. The peak at ~ 13.015 X 10 6 km is ~ 3 times lower for retrograde particles than for prograde ones 
in Model 1. In Model 2 the behaviour is opposite and the rj e value is increased by the same factor of ~ 3. Phoebe's evolution has been 
protected from the effects of these perturbations by its inclination value, which put the satellite into a dynamically safe region. 



10 6 years. It is possible that Narvi and S/2004 S18 had 
similar orbits in the past which later diverged due to the 
inclination-eccentricity link. 

Some of our candidate families merge at higher values of the 
velocity dispersion forming bigger groups we called clusters. 
The most relevant one is that termed as cluster A in table 
[3] It is made of two three-body families and an individual 
satellite and it is defined at a velocity cut-off of 150 m/s. At 
a velocity cut-off of 202 m/s cluster A merges with the two- 
body family related to S/2004 S15 forming cluster B. Con- 
firming these dynamical groups by comparing their colour 
indices is a difficult task because of the limited amount of 
data available in the literature. The only spectrophotomet- 
ry data concerning Saturn's retrograde satellites are those 
of Phoebe and Ymir, which, according to our analysis, are 
separated by a velocity dispersion of more than 600 m/s. 



Phoebe appears to have colours not compatible with any 
other irregular satellite of the system, supporting our claim 
that Phoebe is not related to the rest of Saturn's present 
population of irregular satellites. 

The situation looks more favourable for prograde families: 
the colours of three members of the possible Albiorix family 
are available, with two sets of data for Albiorix itself. The 
colour data are reported in table [4] with the corresponding 
la errors. These data seem to be compatible with the hy- 
pothesis of a common origin of the group at a 3cr level. 



6 CONCLUSIONS 

The aim of this work was to investigate the dynamical 
and collisional nature of the Saturn system of irregular 
satellites. We analysed the secular dynamical evolution of 
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Table 3. Dynamical clustering of Saturn's irregular satellites. The term cluster 
used in the table refers to the merging of previously reported families. 



Family name 



Family members 



Dispersion 



Kiviuq 
Albiorix 
Siarnaq 
Siarnaq + Albiorix 
Prograde 



S/2004 S15 
Mundilfari 
S/2006 S2 
S/2004 S10 
S/2004 S8 
Narvi 
Ymir 

Cluster A 

Cluster B 
Cluster C 

Cluster D 
Retrograde - Phoebe 
Retrograde 



Prograde Satellites 

Kiviuq, Ijiraq 102 m/s 

Albiorix, Erriapo, Tarvos, S/2004 Sll 129 m/s 

Paaliaq, Siarnaq 314 m/s 

Albiorix & Siarnaq families 443 m/s 

All prograde satellites 532 m/s 

Retrograde Satellites 

S/2004 S15, S/2006 SI 114 m/s 

Mundilfari, S/2004 S13, S/2004 S17 116 m/s 

S/2006 S2, S/2006 S3 132 m/s 

S/2004 S10, S/2004 S12, S/2004 S14 144 m/s 

S/2004 S8, S/2006 S5, S/2004 S16 168 m/s 

Narvi, S/2004 S18 200 m/s 

Ymir, S/2006 S2 family, 259 m/s 
S/2004 S7, S/2006 S7 

Mundilfari family, S/2006 S6, 150 m/s 
S/2004 S10 family 

Cluster A, S/2004 S15, S/2006 SI 202 m/s 

Cluster B, S/2004 S8 family, 240 m/s 
S/2004 S9, S/2006 S4 

Cluster C, Ymir family 267 m/s 

All retrograde satellites except Phoebe 315 m/s 

All retrograde satellites 658 m/s 




11.9678 14.9598 17.9517 20.9437 
Semimajor Axis (10 6 km) 



23.9357 



Table 4. Colour indices of three candidate members of 
the Albiorix family: for each colour index the Icr error 
range is reported ( |Grav et al.|2 003). 



Satellite 



V 



V - R 



Tarvos 0.77 ±0.12 0.57 ±0.09 0.88 ±0.11 

Albiorix 0.89 ±0.07 0.50 ± 0.05 0.91 ± 0.05 

0.98 ±0.07 0.47 ±0.04 0.92 ± 0.04 

Erriapo 0.83 ± 0.09 0.49 ± 0.06 0.61 ±0.12 



Figure 27. Mean eccentricities of the retrograde test particles 
computed with Model 2 for different values of inclination. The 
test particles were on initially circular orbits. Most of the irregular 
satellites on retrograde orbits lie in the range of inclination where 
the increase in eccentricity is limited on a 10 6 years timescale. 
Narvi family is a different case. Narvi spends a significant fraction 
of its dynamical evolution in the excited region while S/2004 S18 
is inside it most of the time. As it can be argued from the plot, 
the subsequent divergent evolutions could be at the origin of the 
high velocity dispersion within the group, assuming that the two 
satellites originated from a common parent body. 



the satellites on a timespan of 10 8 years, computing their 
mean orbital elements. We found evidences of resonant 
and chaotic behaviours in the motion of about two third 
of the satellites. We also explored the dynamical features 
of the phase space close to them and verified the influence 
of Titan and Iapetus as well as of the Great Inequality in 
shaping the satellite system. 

In this paper we have also verified that the present orbital 
structure is long-lived against collisions but we found 
indications that in the past a more intense collisional 
activity could have taken place, mainly due to Phoebe's 
sweeping effect. By considering the impact rates of the 
present population, we deduce that the original population 
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could have been at least ~ 30% more abundant. We also 
suggest that the absence of prograde and low inclination 
(~ 170° — 180°) retrograde irregular satellites in the region 
encompassed between 10.47 x 10 6 km and 14.96 x 10 6 
km is a by-product of the sweeping effect of Phoebe. It is 
less clear if the absence of retrograde satellites with lower 
inclinations (i.e. high velocities ejecta from Phoebe or other 
captured bodies) in the same radial region could be due to 
the same reason or if it is a primordial structure related to 
the capture mechanism. 

We also found evidences of dynamical groupings among 
prograde and retrograde satellites, even if their interpreta- 
tion in terms of families is not straightforward due to the 
effects of chaotic and resonant evolution. By applying the 
HCM algorithm and assuming a velocity cut-off of about 
200 m/s, we retrieved two candidate families between the 
prograde and six between the retrograde satellites. Some 
of these families merge in bigger clusters at higher but 
possibly still acceptable velocity cut-offs. This might be an 
indication that the system suffered intense post-breakup 
collisional evolution which could have dispersed the original 
families. 
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